This section is intended to introduce various aspects of the art, which may be associated with exemplary embodiments of the present invention. This discussion is believed to assist in providing a framework to facilitate a better understanding of particular aspects of the present invention. Accordingly, it should be understood that this section should be read in this light, and not necessarily as admissions of prior art.
An important goal of seismic prospecting is to accurately image subsurface structures commonly referred to as reflectors. Seismic prospecting is facilitated by obtaining raw seismic data during performance of a seismic survey. During a seismic survey, seismic energy is generated at ground level by, for example, a controlled explosion, and delivered to the earth. Seismic waves are reflected from underground structures and are received by a number of sensors referred to as geophones. The seismic data received by the geophones is processed in an effort to create an accurate mapping of the underground environment. The processed data is then examined with a goal of identifying geological formations that may contain hydrocarbons.
Seismic energy that is transmitted in a relatively vertical direction into the earth is the most likely to be reflected by reflectors. Such energy provides meaningful information about subsurface structures. However, the seismic energy may be undesirably diffused by anomalies in acoustic impedance that routinely occur in the subsurface environment. Diffusion of seismic energy during a seismic survey may cause subsurface features to be incorrectly represented in the resulting seismic data.
Acoustic impedance is a measure of the ease with which seismic energy travels through a particular portion of the subsurface environment. Those of ordinary skill in the art will appreciate that acoustic impedance may be defined as a product of density and seismic velocity. Acoustic impedance is typically referred to by the symbol Z.
Seismic waves attenuate for a variety of reasons as they travel in a subsurface environment. A quality metric (sometimes referred to a quality factor) Q is typically used to represent attenuation characteristics of underground formations. In general, Q is inversely proportional to seismic signal attenuation and may range from a value of zero to infinity. More specifically, Q is a dimensionless quality factor that is a ratio of the peak energy of a wave to the dissipated energy. As waves travel, they lose energy with distance and time due to spherical divergence and absorption. Such energy loss must be accounted for when restoring seismic amplitudes to perform fluid and lithologic interpretations, such as amplitude versus offset (AVO) analysis. Structures with a relatively high Q value tend to transmit seismic waves with little attenuation. Structures that tend to attenuate seismic energy to a greater degree have lower Q values.
Q values associated with subsurface structures are used to mathematically alter seismic data values to more accurately represent structures in the subsurface environment. This process may be referred to as “Q migration” by those of ordinary skill in the art. During Q migration, a seismic data value representing travel of seismic energy through a subsurface structure having a relatively low Q value may be amplified and broadened in spectrum to a greater degree than a data value representing travel of seismic energy through a subsurface structure having a relatively high Q value. Altering the amplitude and phase of data associated with low Q values takes into account the larger signal attenuation that occurs when seismic energy travels through structures having a relatively low Q value.
FWI is a partial-differential-equation-constrained optimization method which iteratively minimizes a norm of the misfit between measured and computed wavefields. Seismic FWI involves multiple iterations, and a single iteration can involve the following computations: (1) solution of the forward equations, (2) solution of the adjoint equations, and (3) convolutions of these forward and adjoint solutions to yield a gradient of the cost function. Note that for second-order optimization methods, such as Gauss-Newton, the (4) solution of the perturbed forward equations is also required. A more robust mathematical justification for this case can be found, for example, in U.S. Patent Publication 2013/0238246, the entire content of which is hereby incorporated by reference.
A conventional first-order form of the linear visco-acoustic wave equations for simulating waves in attenuating acoustic media is:
                                                                        ∂                p                                            ∂                t                                      +                          κ              ⁢                              ∇                                  ·                  v                                                      +                                          ∑                                  l                  =                  1                                L                            ⁢                                                ϕ                  l                                ⁢                                  m                  l                                                              =                      s            p                          ,                                  ⁢                                                            ∂                v                                            ∂                t                                      +                                          1                ρ                            ⁢                              ∇                p                                              =                      s            v                          ,                                  ⁢                                                            ∂                                  m                  l                                                            ∂                t                                      +                                          κα                l                            ⁢                              ∇                                  ·                  v                                                      +                                          ϕ                l                            ⁢                              m                l                                              =          0                ,                            (        1        )            with appropriate initial and boundary conditions for pressure p, velocity v, and memory variables ml. Note that∇=divergence operator,κ=unrelaxed bulk modulus
      (                            lim                      ω            ->            ∞                          ⁢                  κ          ⁡                      (            ω            )                              ->      κ        )    ,ρ=mass density,v=velocity (v={vx vy vz}T in 3D space),p=pressure,ml=memory variable for mechanism l,sp=pressure source,sv=velocity source,
      ϕ    l    =                    1                  τ                      σ            ⁢                                                  ⁢            l                              ⁢                          ⁢      and      ⁢                          ⁢              α        l              =          (              1        -                              τ                          σ              ⁢                                                          ⁢              l                                            τ                          ϵ              ⁢                                                          ⁢              l                                          )      where relaxation parameters τϵl and τσl may be determined by equation (2) for a given quality factor profile.Note that continuous scalar variables are denoted by italicized characters and vector and matrices are denoted by bold non-italicized characters throughout this document.
                                          Q                          -              1                                ⁡                      (                          x              ,              ω                        )                          =                                            ∑                              l                =                1                            L                        ⁢                                          ω                (                                                                            τ                                              ϵ                        ⁢                                                                                                  ⁢                        l                                                              ⁡                                          (                      x                      )                                                        -                                                            τ                                              σ                        ⁢                                                                                                  ⁢                        l                                                              ⁡                                          (                      x                      )                                                                                                  1                +                                                                            ω                      2                                        ⁡                                          (                                                                        τ                                                      σ                            ⁢                                                                                                                  ⁢                            l                                                                          ⁡                                                  (                          x                          )                                                                    )                                                        2                                                                                        ∑                              l                =                1                            L                        ⁢                                          1                +                                                      ω                    2                                    ⁢                                                            τ                                              ϵ                        ⁢                                                                                                  ⁢                        l                                                              ⁡                                          (                      x                      )                                                        ⁢                                                            τ                                              σ                        ⁢                                                                                                  ⁢                        l                                                              ⁡                                          (                      x                      )                                                                                                  1                +                                                                            ω                      2                                        ⁡                                          (                                                                        τ                                                      σ                            ⁢                                                                                                                  ⁢                            l                                                                          ⁡                                                  (                          x                          )                                                                    )                                                        2                                                                                        (        2        )            
where
Q=quality factor,
τϵl=strain relaxation time of mechanism l in SLS model,
τσl=stress relaxation time of mechanism l in SLS model,
x=spatial coordinate,
ω=frequency,
L=number of relaxation mechanisms used in the SLS model.
Conceptually, the quality factor Q represents the ratio of stored to dissipated energy in a medium. The strain and stress relaxation times are determined to best fit the desired quality factor distribution over the frequency band.
Full wavefield inversion (FWI) methods based on computing gradients of an objective function with respect to the parameters are often efficiently implemented by using adjoint methods, which have been proved to outperform other relevant methods, such as direct sensitivity analyses, finite differences or complex variable methods.
The continuous adjoint of the conventional visco-acoustic system (Equations (1)) is
                                                                        ∂                                  p                  _                                                            ∂                t                                      +                          ∇                              ·                                  (                                                            1                      ρ                                        ⁢                                          v                      _                                                        )                                                              =                                    ∂              ℱ                                      ∂              p                                      ,                                  ⁢                                                            ∂                                  v                  _                                                            ∂                t                                      +                          ∇                              (                                  κ                  ⁢                                      p                    _                                                  )                                      +                                          ∑                                  l                  =                  1                                L                            ⁢                              ∇                                  (                                                            κα                      l                                        ⁢                                                                  m                        _                                            l                                                        )                                                              =                                    ∂              ℱ                                      ∂              v                                      ,                                  ⁢                                                            ∂                                                      m                    _                                    l                                                            ∂                t                                      +                                          ϕ                l                            ⁢                              p                _                                      +                                          ϕ                l                            ⁢                                                m                  _                                l                                              =          0                ,                            (        3        )            
where
p=adjoint pressure,
v=adjoint velocity,
ml=adjoint memory variable for mechanism l, and
∂/∂p and ∂/∂v are derivatives of the objective function with respect to the pressure and velocity respectively.
A common iterative inversion method used in geophysics is cost function optimization. Cost function optimization involves iterative minimization or maximization of the value of a cost function (θ) with respect to the model θ. The cost function, also referred to as the objective function, is a measure of the misfit between the simulated and observed data. The simulations (simulated data) are conducted by first discretizing the physics governing propagation of the source signal in a medium with an appropriate numerical method, such as the finite difference or finite element method, and computing the numerical solutions on a computer using the current geophysical properties model.
The following summarizes a local cost function optimization procedure for FWI: (1) select a starting model; (2) compute a search direction S(θ); and (3) search for an updated model that is a perturbation of the model in the search direction.
The cost function optimization procedure is iterated by using the new updated model as the starting model for finding another search direction, which will then be used to perturb the model in order to better explain the observed data. The process continues until an updated model is found that satisfactorily explains the observed data. Commonly used local cost function optimization methods include gradient search, conjugate gradients, quasi-Newton, Gauss-Newton and Newton's method.
Local cost function optimization of seismic data in the acoustic approximation is a common geophysical inversion task, and is generally illustrative of other types of geophysical inversion. When inverting seismic data in the acoustic approximation, the cost function can be written as:
                                          ℱ            ⁡                          (              θ              )                                =                                    1              2                        ⁢                                          ∑                                  g                  =                  1                                                  N                  g                                            ⁢                                                ∑                                      r                    =                    1                                                        N                    r                                                  ⁢                                                      ∑                                          t                      =                      1                                                              N                      t                                                        ⁢                                      W                    ⁡                                          (                                                                                                    ψ                            calc                                                    ⁡                                                      (                                                          θ                              ,                              r                              ,                              t                              ,                                                              w                                g                                                                                      )                                                                          -                                                                              ψ                            obs                                                    ⁡                                                      (                                                          r                              ,                              t                              ,                                                              w                                g                                                                                      )                                                                                              )                                                                                                          ,                            (        4        )            
where
(θ)=cost function,
θ=vector of N parameters, (θ1, θ2, . . . θN) describing the subsurface model,
g=gather index,
wg=source function for gather g which is a function of spatial coordinates and time, for a point source this is a delta function of the spatial coordinates,
Ng=number of gathers,
r=receiver index within gather,
Nr=number of receivers in a gather,
t=time sample index within a trace,
Nt=number of time samples,
W=norm function (minimization function, e.g. for least squares function (x)=x2),
ψcalc=calculated seismic data from the model θ,
ψobs measured seismic data (pressure, stress, velocities and/or acceleration).
The gathers, data from a number of sensors that share a common geometry, can be any type of gather (common midpoint, common source, common offset, common receiver, etc.) that can be simulated in one run of a seismic forward modeling program. Usually the gathers correspond to a seismic shot, although the shots can be more general than point sources. For point sources, the gather index g corresponds to the location of individual point sources. This generalized source data, ψobs can either be acquired in the field or can be synthesized from data acquired using point sources. The calculated data ψcalc on the other hand can usually be computed directly by using a generalized source function when forward modeling.
FWI attempts to update the discretized model θ such that (θ) is a minimum. This can be accomplished by local cost function optimization which updates the given model θ(k) as follows:θ(i+1)=θ(i)+γ(i)S(θ(i)),  (5)where i is the iteration number, γ is the scalar step size of the model update, and S(θ) is the search direction. For steepest descent, S(θ)=−∇θ(θ), which is the negative of the gradient of the misfit function taken with respect to the model parameters. In this case, the model perturbations, or the values by which the model is updated, are calculated by multiplication of the gradient of the objective function with a step length γ, which must be repeatedly calculated. For second-order optimization techniques, the gradient is scaled by the Hessian (second-order derivatives of objective function with respect to the model parameters). The computation of ∇θ(θ) requires computation of the derivative of (θ) with respect to each of the N model parameters. N is usually very large in geophysical problems (more than one million), and this computation can be extremely time consuming if it has to be performed for each individual model parameter. Fortunately, the adjoint method can be used to efficiently perform this computation for all model parameters at once (Tarantola, 1984).
FWI generates high-resolution property models for prestack depth migration and geological interpretation through iterative inversion of seismic data (Tarantola, 1984; Pratt et al., 1998). With increasing computer resources and recent technical advances, FWI is capable of handling much larger data sets and has gradually become affordable in 3D real data applications. However, in conventional FWI, the data being inverted are often treated as they were collected in an acoustic subsurface medium, which is inconsistent with the fact that the earth is always attenuating. When gas clouds exist in the medium, the quality factor (Q) which controls the attenuation effect plays an important role in seismic wave propagation, leading to distorted phase, dim amplitude and lower frequency. Therefore, conventional acoustic FWI does not compensate the Q effect and cannot recover amplitude and bandwidth loss beneath gas anomalies.
Visco-acoustic FWI, on the other hand, uses both the medium velocity and the Q values in wave field propagation. Thus, the Q-effect is naturally compensated while wave-front proceeds. In some cases, where shallow gas anomalies overlay the reservoir, severe Q-effect screens off signals and causes cycle skipping issue in acoustic FWI implementation. Consequently, a visco-acoustic FWI algorithm and an accurate Q model are highly preferred.
The Q values, however, are not easy to determine. Among many approaches, ray-based refraction or reflection Q tomography has been largely investigated. In field data applications, however, Q tomography is a tedious process and the inversion is heavily depending on how to separate attenuated signals from their un-attenuated counterparts. In recent years, wave-based inversion algorithms such as FWI have been proposed to invert for Q values. Theoretically, such wave-based methods are more accurate. However, the velocity and the Q inversion may converge at a difference pace and there might be severe energy leakage between velocity and Q gradient so that the inversion results are not reliable.
Zhou et al., (2014) describes how to use acoustic FWI for velocity inversion and then how to use the FWI inverted velocity model for Q inversion. However, ray-based Q tomography is time consuming and they did not conduct a real visco-acoustic waveform inversion.
Bai et al., (2014) applied visco-acoustic FWI for velocity inversion, however, they also need to use visco-acoustic FWI to invert for the Q model. As commonly regarded by the industry, such a waveform Q inversion is very unstable. The errors in velocity inversion may easily leak into Q inversion, and vice versa. In addition, this method is hard to do target-oriented Q inversion and Q's resolution and magnitude remain an issue to be solved.